# -*- coding: utf-8 -*-
"""
Created on Sat Sep 17 15:07:22 2016
把手机CBJJAI
@author: zhiqiang
"""

path = r"D:\2016试题\B\B题附件"
fileName = "genotype.dat"
f = open(path+'\\'+fileName)
genData = f.readlines()
j = 1
result = {}
weiDian = genData[0]#表头
weiDian = weiDian.replace('\n','')
weiDian = weiDian.split(' ')
length = len(weiDian)
for gen_i in genData:
    if j==1:
        j+=1
        continue
    gen_i = gen_i.replace('\n','')
    gen_split = gen_i.split(' ')
    for i in range(length):
        result.setdefault(weiDian[i],[])
        result[weiDian[i]].append(gen_split[i])
f.close()
#以上代码会产生一个result字典

#获取性状值
pheno_type_file ='phenotype.txt'
f = open(path+'\\'+pheno_type_file)
pheno_type = f.readlines()
pheno_type_Data=[]
for pt in pheno_type:
    pt = pt.replace('\n','')
    pheno_type_Data.append(pt)


def getRating(result):
    res = {}
    for k in result.keys():
        res.setdefault(k,{'jiankang':{},'bujiankang':{}})
        i = 1
        for vi in result[k]:
            if i<=500:
                res[k]['jiankang'].setdefault(vi,0)
                res[k]['jiankang'][vi] += 1
            else:
                res[k]['bujiankang'].setdefault(vi,0)
                res[k]['bujiankang'][vi] += 1
            i += 1
    return res
    
ress = getRating(result)#运行结果会产生一个{'bujiankang': {'AA': 6, 'GA': 99, 'GG': 395},
 #'jiankang': {'GA': 111, 'GG': 389}}

def getGenWDiffThan_n(ress,n):
    #此函数返回一个和n相差的数组
    reslist = []
    for key in ress.keys():
        if  len(ress[key]['bujiankang'].keys())!=len(ress[key]['jiankang'].keys()):
            reslist.append(key)
            continue
        for k in ress[key]['bujiankang'].keys():
            if abs(ress[key]['bujiankang'][k]-ress[key]['jiankang'][k])>=n and key not in reslist:
                reslist.append(key)
    return reslist
resslist_50 = getGenWDiffThan_n(ress,50)#调用返回resslist_50是相差50的位点


def transformPrefs(prefs):
    '''{'bujiankang': {'CC': 281, 'CT': 179, 'TT': 40},
        'jiankang': {'CC': 333, 'CT': 145, 'TT': 22}}
        转化为：{'CC': {'bujiankang': 281, 'jiankang': 333},
             'CT': {'bujiankang': 179, 'jiankang': 145},
     'TT': {'bujiankang': 40, 'jiankang': 22}}'''
    #此函数转换位点统计量
    result = {}
    for p in prefs:
        for it in prefs[p]:
            result.setdefault(it,{})
            result[it][p] = prefs[p][it]
    return result

first_end=['rs12036216','rs880801','rs4391636','rs7368252','rs7522344','rs2250358','rs707472','rs2273298','rs5746051','rs2999878','rs3013045','rs4646092','rs2143810','rs1883567','rs2807345','rs932372','rs15045','rs7543405','rs9426306','rs12145450']
#['rs3131969'	,'rs28576697','rs6694632','rs2488991','rs12757754','rs11260556','rs11488462','rs2294489','rs1014988','rs7407','rs2748987','rs2460001','rs262687','rs12082516','rs2645065','rs2843160','rs3001344','rs2247308']
def getBatesP(ress,weiDianList=first_end):
    ans = {}
    for keys in weiDianList:#取出特征位点
        tranans = transformPrefs(ress[keys])
        ans.setdefault(keys,tranans)
    return ans

batesP = getBatesP(ress)

def getIndex(test_list,weiDianList=weiDian):
    """将获取到的表头列表，换算为表头下标值"""
    ans = []
    for test in test_list:
        if test not in weiDianList:
            ans.append(-1)
        else:
            ans.append(weiDianList.index(test))
    return ans


def getBatesAns(person,prefs,result):
    p1=1/2
    p2=1/2
    for key in prefs.keys():#获取所有的特征key
        ks = result[key]
        k = ks[person-1]#取出person中在key这个位点上的表征
        p1 = p1*(prefs[key][k]['bujiankang']/500)
        p2 = p2*(prefs[key][k]['jiankang']/500)
    if p1>p2:
        p=1
    else:
        p=0
    return(p)

#'rs12133956','rs3765695','rs6658098','rs7368252','rs2301461','rs10779763','rs590368','rs1888759','rs2807345','rs12752833','rs7543405','rs9426306','rs12036216'
jk_bjk = [getBatesAns(i,batesP,result) for i in range(1,1001)]
jk = jk_bjk[0:500]
bjk = jk_bjk[500:1000]
print('朴素贝叶斯分类预测健康人群准确率',1-sum(jk)/500)
print('朴素贝叶斯分类预测不健康健康人群准确率',sum(bjk)/500)
print('朴素贝叶斯分类预测准确率：',(500-sum(jk)+sum(bjk))/1000)